Surface Hardening and Self-Organized Fractality Through Etching of Random Solids 
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When a finite volume of etching solution is in contact with a disordered solid, complex dynamics 
of the solid-solution interface develop. If the etchant is consumed in the chemical reaction, the 
dynamics stop spontaneously on a self-similar fractal surface. As only the weakest sites are corroded, 
the solid surface gets progressively harder and harder. At the same time it becomes rougher and 
rougher uncovering the critical spatial correlations typical of percolation. From this, the chemical 
process reveals the latent percolation criticality hidden in any random system. Recently, a simple 
minimal model has been introduced by Sapoval et al. to describe this phenomenon. Through 
analytic and numerical study, we obtain a detailed description of the process. The time evolution 
of the solution corroding power and of the distribution of resistance of surface sites is studied in 
detail. This study explains the progressive hardening of the solid surface. Finally, this dynamical 
model appears to belong to the universality class of Gradient Percolation. 



PACS numbers: 64.60Ak, 81.65Cf, 68.35Bs 



I. INTRODUCTION 

Corrosion of solids has major economical consequences 
jl],U . It is also interesting from the point of view of the- 
oretical physics of random systems . 

The comprehension of the basic physical mechanisms 
involved in corrosion implies the study of the dynamical 
evolution of the corrosion process and that of the mor- 
phological features of the corroded surface. 

This paper presents a detailed study of a minimal 
model inspired by recent experiments on pit corrosion 
of aluminum thin films by an appropriate etching solu- 
tion . This two dimensional model is a simplified etch- 
ing model. It was first introduced in where a prelim- 
inary numerical study has been developed. It provides a 
simple description for the action of a finite volume of a 
corroding solution on the surface of a disordered solid. 

When an etching solution is in contact with an initially 
flat surface of a disordered solid, it starts to corrode its 
weakest regions and the surface gets "harder" . However, 
at the same time, new regions are discovered which con- 
tains weak elements. Depending on the corrosion reac- 
tion mechanism, different situations for this hardening 
process can occur. 

Often the corrosive power of the solution is propor- 
tional to an etchant concentration. If the etchant is con- 
sumed in the reaction, then the corrosive power of a finite 
volume of solution decreases during the time evolution 
of the process. As the solid surface gets "harder and 
harder", and the corroding power of the solution gets 
"weaker and weaker" , the corrosion process stops spon- 
taneously in a finite time interval. At this moment all the 
surface sites are "too hard" to be etched by the solution. 



It is this phenomenon which is studied both numeri- 
cally and analytically in this paper. 

A most interesting aspect of this kind of dynamical 
corrosion is that the final surface has a fractal geometry, 
showing that the corrosion mechanism itself uncovers the 
spatial correlations among the strong sites belonging to 
the solid. This is why this phenomenon is intimately re- 
lated to percolation properties of random systems. In 
that sense this kind of corrosion reveals a "latent" criti- 
cality embedded in any random system. 

The model reproduces qualitatively the same phe- 
nomenology observed experimentally || . The dynamical 
evolution can be divided into two different regimes: 

1. In the first (smooth) regime, the corrosion is well di- 
rected and the front becomes progressively rougher 
and rougher. In our model this regime does not 
depend on the details of the discretization chosen, 
not even on the fundamental geometrical features 
of the lattice, like the embedding space dimension 
or the lattice coordination number. 

2. In the second regime, the correlations revealed by 
the hardening process become important: the dy- 
namics becomes locally isotropic generating a frac- 
tal front. This corresponds to a critical regime, 
directly related to the static percolation transition 
on the same lattice. 

The hardness of the final interface, which is related 
to the final corrosion power of the solution, depends on 
the external parameters as the volume of the solution 
itself and the system size. When the volume of the solu- 
tion is not too large, one observes a geometrical scaling 
regime. This regime corresponds to the scaling regime 
of a static percolation model known as "Gradient Perco- 
lation" . When the volume is increased, the correlation 
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FIG. 1. Sketch of the etching dynamics in a square lat- 
tice: the sites 2, 3, 5 are etched at the first time-step as their 
resistances are lower than p(0). Consequently the number 
of etchant particles in the solution decreases by 3 units. At 
t = 1, the new interface sites are 7,8, 10 if the solution can 
etch only the first nearest neighbors solid sites; or the whole 
second layer if the solution can also etch the second nearest 
neighbors in a diagonal direction. 

length grows to reach the system size. Above this limit 
the finite size effects dominate the behavior, and we do 
not study here this case. 

II. THE MODEL 

We first recall the two-dimensional etching model in- 
troduced in ||. Its schematic is shown in Fig. [l]: 

• The solid is represented as a site lattice (triangu- 
lar or square), of linear width L and, eventually 
infinite depth. 

• A random number £ [0, 1] (extracted from the 
flat probability density function 7To(r) = 1 for r G 
[0, 1]) is assigned to each solid site i, representing its 
resistance to the etching by the solution, does 
not depend on time (quenched disorder), and on 
the site environment. 

• The etching solution has a volume V and is ini- 
tially in contact with the solid through the bottom 
boundary (see Fig. [I]) . It contains an initial number 
N et (0) of dissolved etchant molecules. 

Consequently the initial concentration C (0) of etchant 
in the solution is given by: C(0) = N et (0)/V. Calling 
N et (t) the number of etchant molecules at time t, C(t) — 
N e t(t)/V. At each time-step, the "etching power" of the 
solution (i.e. the average "force" exerted by the solution 
on a solid surface particle) is supposed to be proportional 
to C(t) : p(t) = TC(t). Hereafter the assumption T = 1 
is made, without loss of generality. It implies C(t) = p(t). 
At time-step t, all the interface sites with < p(t) are 
dissolved and a particle of etchant is consumed for each 
of these corroded solid sites. 

Let us call n(t) the number of dissolved solid sites 
at time t. One can express many important dynamical 
quantities through n(t), or its time-integral N(t), that is 




FIG. 2. Typical process represented at two intermediate 
time steps, and at the final one. The solid is represented in 
grey, the solution in white, and the finite size solid clusters 
detached by the solution in black. The final solid surface is 
found to be fractal up to a characteristic scale a. 

the total number of corroded solid sites up to time t. The 
number of etchant particles in the liquid will decrease as: 

N et (t + 1) = N et {t) - n(t) = N et (0) - N(t) , (1) 

and consequently the etching power of the solution is: 

p(t + l)=p { t)-^f=p(0)-^. (2) 

Note that, as p(t + 1) < p(t), a site having resisted to 
etching at a certain time-step will resist forever. Conse- 
quently, the part of the solid surface which can be etched 
at time-step t + 1 is restricted to the sites which have 
been just uncovered by the etching process at time t. Wc 
call this subset of surface the "active" part of the sur- 
face. After a given time-step, all the solid sites which 
have been previously explored by the solution are defi- 
nitely "passive" . However it may happen that "passive" 
sites are disconnected from the bulk at a later time-step 
if they are connected to the solid by weak sites. 

A. Phenomenological description of the dynamics 

A typical process at two intermediate times, and at 
the final time-step, is shown in Fig. ^[ Some finite solid 
clusters are detached from the "infinite" solid by the cor- 
rosion process. Consequently, at any time, the "global 
surface" of the system is composed by both the finite 
clusters surfaces and the surface of the infinite solid , 
which will be called the "corrosion front". Note that, in 
order to have a meaningful geometrical and physical defi- 
nition of the solution space and of the connected solid re- 
gions (and then of the corrosion front), one has to use the 
so-called "dual" connection rules for solution and solid, 
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respectively [Q. For example, on the square lattice, if 
the solution etches both first and second nearest neigh- 
bors, only first nearest neighbor solid sites should be con- 
sidered as connected. On the other side, if the solution 
etches only first nearest neighbor, both first and second 
nearest neighbors solid sites should be considered as con- 
nected. On the triangular lattice, if the solution etches 
first nearest neighbors, liquid and solid sites are consid- 
ered to be connected both by first nearest neighbors only. 
Two remarks should be made: 

• The corrosion front keeps quite smooth at the be- 
ginning of the dynamics (first snapshot of Fig. |^ ). 
It becomes very irregular only towards the end of 
the dynamics (third snapshot of Fig. |J) when p(t) 
is close to the percolation threshold p c on the same 
lattice ||. 

• The "active" part of the global surface is essentially 
restricted to the etching front, since a site having 
resisted to the corrosion at a certain time-step will 
resist forever. 

These observations are useful for a first analysis of the 
dynamics. Roughly speaking, if the front advances lin- 
early, the number of solid sites discovered at each time- 
step is L (the number of site in each layer). Hence, 
in this approximation, the number of etched sites is 
n(t) — L p(t). Using this approximation one gets (from 
Eq.l): 

Pit) = p(* " 1) (l " = P(0) (l - |) ' ■ (3) 

This simple prediction is compared with the actual simu- 
lation behavior of p(t) in Fig. ||. The agreement between 
the simple prediction [}] and the initial decay of p(t) is 
very good for values p(t) > p Cl i.e. in the smooth regime 
of the dynamics. When p(t) is close to p c this approxima- 
tion is no more valid and the dynamics enters the critical 
regime. 

A better derivation of Eq. |3[ and a more precise def- 
inition of the two regimes will be given below providing 
a deeper insight on the critical regime of the dynamics, 
when the surface becomes fractal and the dynamics slow 
down and stop. 

Note that the main hypothesis for the derivation of 
Eq. [3] consists in assuming that at each time-step the 
number of new sites checked for corrosion is always L, i.e. 
the whole next solid layer. This is possible if the etching 
does not leave large connected segments of uncorroded 
sites. In fact it is easy to show that the non-etched sites, 
the number of which is approximatively (1 — p(t))L, are 
almost isolated, the average size of a segment of "sur- 
vived" sites being (I) = l/p{t). 

Interestingly the present phenomenological approach 
suggests an analogy between our dynamical etching 
model and a static percolation model known as Gradi- 
ent Percolation p|JlO||. This will be discussed next. 
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FIG. 3. Decay of the corroding power p(t) in a triangular 
lattice with p(0) = 1. The numerical evolution of p(t) (solid 
line) is compared with the simple phenomenologic derivation 
of Eq. |^ (dashed line). When the two curves separate the 
dynamics enters the "Critical Regime" dominated by perco- 
lation effects. 

B. Analogy with Gradient Percolation (GP) 

The Gradient Percolation (GP) problem fl(],[l3| can 
be formulated in the following way: a random number 
Ti G [0, 1] is assigned to each site of a lattice of rr-size L 
and y-size h. A constant gradient of occupation proba- 
bility in the y-direction is then imposed on the lattice: 
p(y) = 1 — y/h = I — y Vp. The occupation rule is 
that in each column the sites is occupied if and only if 
i'i < p(y) ( see Fig. |]). In the first column (y = 0) the oc- 
cupation probability is one, while in the last one (y = h) 
it is zero. These two special layers individuate two per- 
colating clusters in the x direction of occupied (grey) 
and empty (white) sites (Fig. The external frontier 
of the connected occupied cluster is called the gradient 
percolation front fllPf . This front is centered around the 
layer with p(y) equal to the critical percolation threshold 
p c characteristic of the lattice type. The front is frac- 

(GP) 

tal with a dimension D ^ ~ 1.75 up to a finite length 
(front width o~gp) which is a power law of the local gra- 
dient Vp = l/h: 

°gp ~ [Vp]- Q ^ P) (4) 

where af P) ~ 0.57. Note that D { f GP) ~ 7/4 and 
a„ ~ l/Dj P \ For this reason it was assumed that 
Df P is equal to the fractal dimension of the hull of the 
incipient infinite percolating cluster D'j = 7/4 in percola- 
tion theory Jl7],[l8| . The demonstration of the identity of 

(CP) (CP) 

the equivalent of a a = 1 /D \ in percolation theory 
is given in |l7| |. 

In addition, the occupation probabilities of the front 
range in an interval p{y) — p c ± Ap, where Ap scales 
with the gradient as 
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FIG. 4. The Gradient Percolation Model. The numbers r 
are chosen randomly between and 1. Each column has an 
occupation probability p(y) ranging from 1 at the left side to 
on the right side. A site i is occupied if ri < p(y). The 
occupied and empty sites are represented respectively in grey 
and white. Apart from isolated islands and lakes, grey and 
white sites form two distinct connected regions. The bold line 
represents the separation between these two regions. 



Ap ~ [Vp\ 



(5) 



The exponent a^ P ^ is related to a^ P ^ , as Ap ~ a Vp, 
which implies, from Eq. ^, ap GP ' = 1 — a^ P ^ . 

Because of its characteristic properties, GP has pro- 
vided a powerful method to compute percolation thresh- 
oldp c 0. 

In this model, one can associate for each corroded site 
(x, y) the value p(x,y) of the solution etching power at 
the time of corrosion of that site ||] . In this way, a posi- 
tion dependent "field" of occupation probabilities (by the 
solution) is spontaneously generated. This is the physi- 
cal link with GP. In the smooth time regime the succes- 
sive active zones are consecutive solid layers containing 
about L sites. Consequently, p(x, y) depends only on y 

(p(x,y) =p(y))- 

The "active" zone at time t is then the whole layer at 
depth y = t From Eq. one can then write: 



P(y)=p(0)[l 



(6) 



This equation defines a Self-Organized Gradient Perco- 
lation, where the value of the gradient depends on the 
parameter L and V as: 



(Vp). 



L 
V 



(7) 



III. SIMULATIONS AND NUMERICAL RESULTS 

Extensive simulations have been performed, consid- 
ering triangular and square lattices, with first nearest 



neighbour (f.n.n.) and second nearest neighbors (s.n.n.) 
(diagonal) connections for the corrosion process. All sim- 
ulations start with po = 1 > p c in order to observe clearly 
the transition towards a critical regime, when p(t) ~ p c . 
Once po is fixed, the parameter measuring the initial cor- 
roding "force" of the solution is V = N et (0)/po- The 
other fundamental parameter is the transversal size of 
the solid L. All the data presented below refer to 1000 
different realizations of the quenched disorder, for each 
choice of the parameters L and V. 



A. Correlation length and "Phase" diagram 

In order to quantify the statistical properties developed 
by the dynamical process, the average thickness of the 
final corrosion front is measured. If {yi} are the depths 
of the points i belonging to the corrosion front at time t, 
its average thickness can be defined as: 



\ 




where / is the length of the corrosion front. 

The behavior of the final value a at time tf as a func- 
tion of the "natural" gradient L/V is shown in Fig. ^ 
(bottom) for different fixed values of L. Several observa- 
tions can be made: 

• First, for sufficiently large values of L/V (right side 
of Fig. |sj> a follows the scaling behavior 



(L/vy 



(8) 



with a a = 0.57 ± 0.02. This confirms the idea that 
the final features of our dynamical etching model, 
at least in this range of L/V, are in the same univer- 
sality class of GP once the identification L/V ~ Vp 
is done, i.e. a„ 



(GP) 
- Ota 



• Decreasing L/V, cr increases following the previ- 
ous scaling behavior (Eq. ||) until reaching values 
of L/V for which a ~ L. For even smaller values 
of L/V, a deviation from the aforementioned scal- 
ing law is observed. This deviation is characterized 
by a cross-over to a region dominated by boundary 
effects. In this regime a seems to decrease slowly 
together with the gradient L/V, instead of increas- 
ing. 

• Consequently, for a fixed value of L, one can dis- 
tinguish a "strong gradient" process, i.e. for values 
of L/V in which Eq. || holds, and a "weak gra- 
dient" process for values of L/V smaller than the 
cross-over value. The cross-over between the two 
behaviors is marked by a marginal value of L/V 
for which a ~ L. Note that for this value of L/V 
the spatial correlations extends all over the sample. 
Then this is a kind of "critical" value of L/V. 
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FIG. 5. Determination of the different corrosion behaviors. 
Bottom: behavior of a as a function olL/V for several sample 
sizes (several values of L). Top: the behavior of (p/) for the 
smaller sample (L = 500). Note that the maximum of a 
corresponds to the change of sign of (pf) — p c (vertical arrow). 
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FIG. 6. "Phase" diagram in the parameters space (see 
text). 



Moreover one observes that (pf) < p c in the strong gra- 
dient regime and (pf) > p c in the weak gradient regime, 
where (...) means an average over different realizations of 
the disorder with fixed parameters L, V. In this way the 
equality (pf) = p c can be used to identify the marginal 
("critical") value of L/V for a fixed value of L. 

In the upper diagram of Fig. ^], the transition between 
the two regimes for L = 500 is shown. This transition 
corresponds to a value of L/V ~ 2 • 10~ 5 (marked by the 
double arrow crossing the two plots). 

This behavior of a allows to sketch a kind of "phase" 
diagram for our model in the (L, V) parameter space 
(Fig. ||). The "critical" line a ~ L separates the two 
"phases" . Here we use the terminology of phase transi- 
tions because in GP the correlation length is equal to the 
front width a. 

Since in the strong gradient "phase" Eq. || holds, the 
scaling relation for the marginal line is: 



-JXS 



FIG. 7. Typical final corrosion front for a simulation in 
the strong gradient phase. Note that a <C L. Different 
non-overlapping almost independent regions are identified by 
numbers. The number v = L/a of almost independent regions 
is relevant for the study of extremal quantities [IS]. 
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L/V 

FIG. 8. Scaling behavior of p c — {pf ) . Note that, identifying 
L/V with Vp of Gradient Percolation, one obtains the same 
values of the exponent a p — 0.46 ± 0.02. 



L - V— 



(9) 



The relevance of this relation with respect to the ex- 
tensivity of spatial correlations in the "thermodynamic 
limit" is discussed in Appendix A. In the following, we 
deal only with the strong gradient regime, leaving the 
detailed analysis of the weak gradient regime to further 
work. 



B. Strong gradient etching 

In order to study this regime, we focus on simulations 
of sizes L — 3000 and L — 5000, with (pA < p c . Such 
values of L are large enough, and at the same time they 
permit to collect large statistics. A typical corrosion front 
is presented in Fig. |7| where the condition L 3> a is em- 
phasized. Note that on scales larger than a, the corrosion 
front is almost flat. This indicates the statistical inde- 
pendence among non-overlapping regions of the surface 
of linear size larger than a. 

As mentioned earlier, a is described by Eq. ||. Similarly 
to cr, other important properties follow simple scaling 
relations with the gradient L/V ||. 

The distance of the average value (pf) from p c follows 
the scaling law: 
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FIG. 9. Scaling behavior of the average number of sites per 
column in the final corrosion front. 



Pc-(Pf)~{y) , with a p = 0.46 ±0.02 (10) 

as shown in Fig. [|. f] 

Moreover, the average number of corrosion front sites 
per column (L(tf)) /L is found to follow a power law of 
the form: 



percolation This effect can explain the reduced 

fractal dimension (4/3) measured in the real corrosion 
experiments fjj , due to insufficient image resolution. For 
example, on the triangular lattice (where the solution 
etches only f.n.n.), if the resolution does not distinguish 
first and second nearest neighbors, the measured fractal 
dimension is 4/3. 

The average critical time t c , defined by p(t c ) = p Cl 
and the difference between the arrest time tf of the dy- 
namics and t c itself, are measured for different values of 
the gradient L/V . For the first one, the following simple 
behavior is found (see Fig. |TT| ) : 

(t c ) ~ {L/V)~ at " with a tc = 0.998 ± 0.001 . (12) 

As we shall see in the following, this is a direct conse- 
quence of the linear properties of the smooth dynamical 
regime (Eq. ||). Finally, for tf one has (see Fig. [l2]) 



(t f ) - (t c ) ~ {L/Vy at f with a tf « 0.55 . 



(13) 



However, for (tf) — (t c ), a further dependence on L is ob- 
tained (see the inset of Fig. [l2|). In particular, changing 
L with L/V fixed, the quantity (tf) — (t c ) is found to 
depend linearly on In L. This behavior is connected to 
the "extremal" nature of tf and is not studied here |L9| . 



v 



with ai = 0.45 ± 0.02 . (11) 



The fractal dimension Df of the corrosion front was 
measured (up to the scale a) using the box-counting |l2] ] 
algorithm. In this way Df = 1.753 ± 0.005 is found (see 
Fig. [n^). Note that it is compatible with the value 7/4 
of GP. 

In the early papers about this etching process Q Df ~ 
1.62 was measured. This different value was due to finite 
size effects. The present simulations are almost 400 times 
larger than those of those previously reported in || (the 
largest value of the parameters being L — 3 ■ 10 4 and V = 
5 ■ 10 9 ). This achievement is important to assert that the 
exponents characterizing the final corrosion front belongs 
to the universality class of Gradient Percolation. Whilst 
(pf) depends on the lattice geometry, as p c changes, the 
values of the exponents remain the same. 

Nevertheless, note that the measured fractal dimension 
of the corrosion front can be reduced to 4/3, if one does 
not use the right "dual" connectivity criterion introduced 
above. This is the so-called Grossmann-Aharony effect in 



* Appendix A discusses the possibility of obtaining a value 
of (j)f) arbitrarily near p c , remaining in the "strong gradient" 
region of Fig. ^| Starting from a couple (Lo, Vb) in the strong 
gradient phase, one obtain it, for instance, performing the 
limit V -> oo on any line (L/L ) = (V/V ) a with a a /(l + 
a a ) < a < 1. 



C. Scaling relations 

The exponents a a , a p , aj, and Df are not indepen- 
dent. At first, note that, within the present numerical 
precision, 



1 

Df~ 



(14) 



as in GP gg. 

Identifying the width a with the horizontal correlation 
length, the average number within a correlated region 
scales as a Df because of the fractality on smaller scales. 
Since the horizontal size of the solid is L, the average 
number of distinct correlated regions will be L/a. Con- 
sequently, one can write: 



(I(t f ))~-* D f, 
a 



which implies 

</(*/)) 



Df-l 



-a a (D t -\) 



(15) 



From Eqs. |TT|, |lj and one then obtains the following 
scaling relation: 



ai = a a (Df - 1) 



D i 



Df 



(16) 



which is consistent with the measurement of ai in Eq. [ll]. 
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Box Counting 




FIG. 10. Box-counting determination of the fractal dimen- 
sion Df of the corrosion front. The value of Df = 1.753±0.005 
is found fitting for values of e ranging from a few lattice dis- 
tances to the front width a (in this case sigma w 3000, i.e. 
ln(o-) « 8). 
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FIG. 11. Scaling behavior of the critical time t c for which 
p(tc) =pc- (tc) ~ (L/V)~ at ' with a tc = 0.998 ±0.001. 
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FIG. 12. Scaling behavior of (tf) — (t c ) for different values 
of L. Inset: the dependence on L for L/V fixed is presented. 



Exploiting the analogy between L/V and the gradient 
Vp in GP, another interesting relation among exponents 
can be derived. From the relation Ap = p c —pf ~ Vp • a, 
one gets: 



a p = 1 



Df 



(17) 



Note that this implies a p = aj in d — 2. In fact the as- 
sumption that the number of different correlated regions 
scales as L/a is valid only in d = 2. 



IV. DYNAMICAL EQUATIONS AND 
THEORETICAL RESULTS. 



In this section we present an analytical derivation of 
the dynamical evolution of p(t) and the distribution of 
the surface resistances. This time-dependent distribution 
characterizes the evolution of "hardening" properties of 
the surface. To this aim, the histogram h(r, t) is intro- 
duced. The quantity h(r,t)dr measures the number of 
global surface sites with random resistance in the inter- 
val [r, r + dr] at time t. By definition, the number of sites 
in the global surface G(t) is simply the total integral of 
the histogram: 



G{t) 



h(r, t)dr . 



(18) 



On the other side, the number of surface sites being cor- 
roded at time t by the solution will be: 



n(t) 



p(t) 



h(r, t)dr , 



(19) 



as n(t) is the number of sites in the global surface with 
r < p{t). Note that Eq. |l| links h(r,t) directly to p(t) 
through Eq. 0, which can then be rewritten: 



I P(t) Hr,t)dr 



p(t + l)=p{t) 



V 



(20) 



Let us call m(t) the number of active sites at time t + 
1: i.e. the new sites entering the global surface as a 
consequence of the corrosion of the set of n{t) sites. Then 
the set m(t) is the active zone at time t + 1. One can 
define w(t) = m(t)/n(t). Therefore u>(t) is the number 
of new active sites per etched site at time t. As shown 
below, the quantity u)(t) is the fundamental parameter 
relating the "geometry" to the "chemistry" of the system 
at time t. At each time-step one can write 

G(t+l) = G(t)-n(t)+m(t), (21) 

or, using both Eq. [l9| and the definition of u(t): 

rv(t) 

G(t + 1) = G(t) + (u(t) - 1) / h(r,t)dr. (22) 

Jo 
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Considering only sites in [r, r + dr], one can write: 



A. Smooth regime 



h(r, t + 1) = h(r, t) - h(r, t)8(p(t) - r) 
rp(t) 

+u(t) / h{r',t)dr' , 



(23) 



where 8(x) is the Heavyside step-function. In Eq. ^3[ the 
second term in the right-hand side represents the num- 
ber of sites etched at time t (a surface site i is etched 
with probability 1 if rj < p(t) ). The third term is the 
contribution to h(r, t) due to the new active zone. It is 
based on the fact that each new site has completely ran- 
dom resistance to etching; the probability that it belongs 
to the interval [r, r + dr] is simply dr (as 7i"o(r) = 1). 
In principle, knowing the behavior of u>(t), one can solve 
the system given by Eqs. [2^ and 23, characterizing in this 
way the dynamical evolution of the corrosion power and 
of the resistance of the solid surface. 

Before going on with calculations, it is important to 
observe that h(r,0) = L Vr 6 [0,1] (as n {r) = 1 and 
the initial surface is a layer of sites of length L). On the 
other hand for r < p(t), Eq. reduces to 



P(t) 

h{r,t + \) =w(t) I h(r',t)dr' 



(24) 



Eq. and the initial condition h(r, 0) = L imply that 
at each time for r < p(t — 1) h(r,t) is independent on r 
and can be written as: 

t-i 

h(r,t) = L Y[(cj(t r )p(t'j) far r<p(i-l) (25) 

t'=0 

Using this expression in Eq. ^o], the following equation is 
obtained: 



p(t + l)=p(t) 



V 



t'=0 



(26) 



Eq. |26| makes evident the strong dynamical link between 
the geometry (w(t)) and the chemistry (p(t)) of the sys- 
tem. 

In order to examine the calculations further, it is nec- 
essary to make some hypothesis on the behavior of u)(t). 

As previously mentioned, the dynamical evolutions can 
be divided into two regimes: 

1. a first smooth regime, which is referred to the time 
scale at which p(t) is larger than p c ; 

2. a second critical regime, which is referred to the 
time scale at which pit) ~ p c . 

This partition of the dynamics into two regimes is di- 
rectly connected to percolation theory Jl8| , as shown be- 
low 



If one considers all the lattice sites with r < pit) for 
pit) > p c , they form both a set of a few finite size clusters 
and an infinite percolating and homogeneous (not fractal) 
cluster fl8|] . Consequently, the intersection between the 
global solid-solution surface and this set is made of a 
large number of sites. The larger pit), the larger the 
intersection. This intersection is nothing else but the set 
of n(t) sites to be dissolved at that time-step. 

Since n{t) 3> 1 (and then m(t — 1) > nit) ^> 1 also), 
one can use the law of large numbers to relate n(t) to 
m(t - 1): 



i(t) =p(t)m(t-l). 



(27) 



For the same reason one expects small fluctuations 
around these values. From Eq. ^ and the definition of 
one obtains 



m(t) = Luit)pit)mit - 1) . 



(28) 



Because of the percolation properties for p > p c , which 
are related to the previous argument, one expects 



m(t) — m{t — 1) 



m(t - 1) 



< 1 



Hence the relation: 



w(t) 



1 

W) 



(29) 



This equation introduces an important relationship be- 
tween the fundamental "chemical" parameter pit) and 
the "geometrical" and "dynamical" parameter w(t). In- 
serting Eq. ^ in Eq. one recovers Eq. |^, which can 
be rewritten as: 



pit) = p(0) exp - 



t 



(30) 



withr= -l/ln(l-L/V) ~V/L (L>V). FromEq.|^ 
one has t c — V/Llnipo/p c ) ~ (L/V)^ 1 - This confirms 
the numerical result found in Fig. |ll| and expressed by 
Eq. Note that the behavior given by Eq. ^0| is inde- 
pendent on the space dimension and on the coordination 
number of the lattice. This behavior is however valid only 
up to the time at which pit) ~ p c . After this time the 
hypothesis to deduce Eq. |30| and in particular the possi- 
bility of using the law of large numbers, breaks because 
of the geometrical constraints given by the percolation 
properties of random numbers on a lattice. 

Using Eq. [50], one can derive rigorously the shape of 
h(r, t) at any time-step t or of its normalized version 
4>ir,t) (i.e. J Q dr(j>{r,t) = 1). <fiir,t) is obtained by divid- 
ing h(r, t) by Git). Technical calculations are reported 
in Appendix B. We provide here directly the result: 
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(r, t) ~ x 



j for r < p(i) 

x { 1- ^(lnp(O) -Inr) for p{t) < r < p(0) 
1 for r > p(0) 



(31) 



where 



[MtT 1 = 1 + (p(0) - p(*))/ln(p(t)/p(0)) 
= l-^p(0)(l-exp(-t/r)). 



(32) 



Eq. 31 can be rewritten, in term of instead of t, as 
follows: 



0(r,t)~^i(f) x 

for r < p(i) 
Sft farp(*)<r<p(0) (33) 



for r > p(0) 



B. Critical regime 

For p(t) — p c , one has (in the limit L — > oo) a marginal 
critical case in which the set of lattice sites with r < 
p(t) form finite size clusters of any size and an infinite 
fractal percolating cluster. Finally for p(t) < p c , the set 
of lattice sites with r < p(t) forms only finite size clusters. 
For this reason, even if at such time t the intersection 
between the global solid surface and the set of lattice 
sites with r < p{t) is not empty, it becomes depleted 
after a finite number of time-steps. The average number 
of time-steps after which the dynamics stops will be a 
function of the system parameters L and V (Eq. 13). 
At this time the corrosion dynamics stops because the 
intersection between the global surface and the set of 
lattice sites with r < p{t) is empty. This explains why the 
final corrosion front is fractal with a fractal dimension D f 
and a characteristic size (thickness) a. Df = 7/4 is the 
hull fractal dimension of the finite clusters formed by the 
lattice sites with r < p(tf) and a is the characteristic size 
of these clusters. Finally, the same argument explains 
why each exponent, characterizing the above introduced 
scaling relations (apart from those about t c and tf), is 
directly connected to the exponents of GP. 

From the above argument, it is important to note that, 
if not empty, the active zone at a time t > t c is com- 
posed by a small and fluctuating number of sites m(t). 
This implies that also n(t) is small and strongly fluctu- 
ating. Consequently, the arguments developed in dealing 
with the smooth time-regime, based on the law of large 
numbers and small fluctuations, are no longer valid: uj(t) 
becomes a strongly fluctuating quantity. These critical 
fluctuations are related to the fractal morphology of the 
critical phase of percolation. We can say that the ar- 
rest of the etching dynamics is due to one of these big 
fluctuation of lo in which no site of the active zone has 
r < p(t). 



All these features are shown by Fig. |l3|, where uj(t) and 
\/p{t) are shown as functions of time. It is important to 
note that, whereas co(t) is a strong fluctuating quantity 
in the critical time-regime, p{t) is always smooth. In fact 
p(t) can be written (Eq. |^) as 



p(t) = p(0) 



t-i 



fc=0 



Consequently, p(t) can be seen, apart from prefactors, as 
the time integral of n(t), which is a limited function of 
time and then p(t) is continuous. Moreover, Fig. [l3| shows 
that in the critical time-regime the equality uj(t) ~ \/p{t) 
is valid only "in average" : 



1 



(34) 



where uj(t) means the average of u) over a sufficiently 
large time interval around t in the critical time regime. 
In order to justify the smooth behavior of p{t) in spite the 
fluctuations of uj(t), one can use Eq. |2^. In the continuous 
time limit, it can be rewritten as: 



dp(t) L 

— =__p(t)exp 



dt'{\np(t') + \nuj{t')) 



(35) 



Since the fluctuating tu(t) appears in the integral term 
of the exponential of the right-hand side of Eq. |3^, one 
expects p(t) to be regular. 

In order to understand why one has u>(t) ~ l/p{t) also 
in the critical regime, one has to analyze the behavior of 
time average 



t+At/2 
t'=t-At/2 

Using the definition of to(t), one can write: 

t+At/2 

At 



t'=t-At/2 

From Eq. one has directly 



m(t' - 1) m(t') 



n{t') 'm(t'-l) 



(36) 



t+At/2 

J2 [m(t'-l)/n(t')} = l/p(t). 

t'=t-At/2' 

From this observation and supposing that the quantity 
m(t')/m(t' — 1) oscillates symmetrically around 1, one 
expects Eq. [34|. Before proceeding, it is worth to note 
that the variation of p(t) during the critical time-regime, 
as shown by Eq. |To|, is very small. From this one has 
approximatively: 



L0{t) ~ l/p c 



(37) 



Eqs. |34| and 37 are very important because they pro- 
vide both a "physical" and a "geometrical" meaning to 
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1200 



FIG. 13. Comparison between u(t) (dots) and l/p(t) (thick 
line) as functions of time. 



the critical threshold of percolation in a given lattice (an 
analogous relation for Invasion Percolation was found in 
©@). 

In order to clarify the geometrical effect on u>{t), it is 
important to observe that, following Eq. ^9] and Eq. |3(], 
ui(t) would increase to infinity. Clearly this is not pos- 
sible in a finite dimensional lattice with a finite coordi- 
nation number. For instance in a 2d site square lattice 
with f.n.n. connection u) must be always be smaller than 
3. Moreover, as seen above, the percolation theory in- 
troduce a stronger constraint forbidding p(t) to go well 
below p c . ^ Consequently, in order to use Eq. ^6] to 
predict the behavior of p(t), one has to take in account 
both the behavior given by Eq. [m] and these geometri- 
cal constraints. In order to show this, we have made the 
approximation u>(t) — l/p s i m (t) m Eq. |2^, where p s imif) 
is the simulation outcome for p(t) (and then it includes 
automatically the geometrical constraint). The solution 
p(t) obtained by solving numerically Eq. ^6] is very near 
to Psim(t) itself. 

Because of the strong fluctuations, the purely analyti- 
cal study of this critical regime is difficult. For this reason 
we developed only an approximated mean field approach 
by imposing only the geometrical constraint with the fol- 
lowing simple approximation: 



u{t) = 1/pc 



(38) 



in this critical time-regime. This is a kind of mean field 
approximation as the fluctuations of w(t) are neglected. 
Inserting the relation in Eq. E6l one can write: 



^Note that p c is always larger than the inverse of the coordi- 
nation number of the lattice |H| . They are equal only in the 
Bethe lattice without loops. 



p(t+l)=p{t) 



i-£pe-*n*o 



t-1 



V 



t'=0 



(39) 



with the initial condition p{t c ) = p c and t > t c . In order 
to solve Eq. |3^, one can consider the continuous limit 
which is equivalent to take u)(t) — l/p c in Eq. [sH]: 



dp(t) 
dt 



L , \ 
-— p(i)exp 



Pc 



(40) 



This equation can be solved exactly, and this solution is 
well approximated by: 



Pc - p(t) ^ Pc 



2L ex P (\/^(* 



-t P . 



- 1 



CXp (^F(f - t c )) + 1 



(41) 



where t c is the time at which p(t) — p c . From Eq. [Il] we 
can sec that in the t — > oo limit one has 



Pc - Pj 



(42) 



In spite of the rough approximation, we obtain a good 
approximation of the numerical result of Eq. ^ (see also 
Fig. ||), then this "mean field" result provides a good 
approximation. 

Moreover, since the time constant of Eq. fll| is t' ~ 

Y, it is argued that the average time tf at which the 
dynamics spontaneously stop, obeys the following scaling 
law: 



tf - t c 



Also, this is a reasonable approximation with respect to 
the simulation behavior (Eq. |l3| and Fig. |l^) . We can see 
that Eq. [l|] can be interpreted as the product of the Ap in 
single time steps ( Ap ~ L/V) multiplied by the average 
number of time-steps necessary to stop the dynamics tf — 

tr r 



V. FINAL HISTOGRAMS AND SURFACE 
HARDENING 

The surface hardening can be described by considering 
both the histogram of the global surface and that of the 
corrosion front. One first observes that the number of 
corroded sites in the critical regime is much smaller than 
the number of sites which have been etched in the smooth 
regime. This means that the global surface is dominated 
by sites belonging to finite clusters. As a consequence, 
4>{r,tf) is well approximated by the linear regime be- 
havior given in Eq. B7 with t = tf and p(t) — pf (see 
Fig. [lj]). One observes good agreement between theory 
and numerics. 
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FIG. 14. Normalized histogram of site resistances belong- 
ing to the final global surface. Circles represent the result of 
numerical simulations and the solid line the theoretical pre- 
diction detailed in Appendix B. 
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FIG. 15. Normalized histogram of the corrosion front resis- 
tances. One can notice that simple step functions are found 
for any three different lattices studied. The discontinuities 
occur very near the percolation thresholds. 



The global histogram describes the hardening phe- 
nomenon of the global surface which includes finite size 
clusters detached by the etching process at various time- 
steps of the dynamics. The increasing behavior oi(/)(r,tf) 
for r > p c is due to the fact that the majority of surface 
sites belonging to finite islands have been discovered at 
intermediate time-steps when p(t) > p c . This implies 
that their resistance is well above p c . 

The hardening of the corrosion front is described by 
the normalized distribution <Pf(t) of it site resistances. 
It has been measured by numerical simulation for several 
types of lattices. The numerical results are shown in 
Fig-d 

As discussed above, the final front is more resistant to 
etching than the original native surface. This is shown in 
Fig. Fta, giving the histogram of the front resistances. In 



first approximation it is a step function around r = p c . 
This confirms the hypothesis, derived by GP, that the 
final corrosion front corresponds to the hull of percolation 
clusters with p = p$ ■ 

This effect could possibly be used practically in waste 
management problems. Consider for instance a sys- 
tem containing dangerous compounds with random re- 
sistances to etchants present in the environment. If in 
natural circumstances it comes in contact with etchants, 
even with a weak etching power po < p c , then dangerous 
materials can be diffused in the environment. 

But one can think to apply to the system a previous 
etching treatment with p > p c , in an artificially con- 
trolled situation. In this case the final surface contains 
only sites with r > p c . Then the treated surface will 
resist forever to any further natural attack with p < p Cl 
with no danger of leaks in nature. This could be called a 
random "Darwinian" selection of a strong surface : once 
selected with a finite solution with p > p c the surface 
resists for ever to further etching with pq < p c whatever 
the volume of the solution. 



VI. DISCUSSION 

Several properties and limitations of this model should 
be discussed. This depends on the nature of what has 
been described here as a "site" . One can think of a site 
as being an atom or a small group of atoms. For example, 
if one considers a random solid (like a glass) the different 
local random environments will cause random local rates 
for atomic dissolution rather than random probabilities. 
The possible application of the above model is then re- 
stricted to situations where the choice of suitable time 
intervals makes it possible to separate "very resistant" 
and "very weak" sites. This implies a transposition of 
the distribution of local rates into a distribution of sites 
resistances. 

But a site could also be a semi-macroscopic entity like 
a small crystallite protected by a randomly resistive sur- 
face. In the case of corrosion experiments by Balazs || 
it is believed that randomness may be attributed to the 
random nature of the oxide layer which spontaneously 
grows on previously etched aluminum crystallites. The 
disorder studied herein occurs if the random resistances 
to etching appear just after oxidation of newly discovered 
crystallites. Although the disorder appears dynamically, 
once created it is quenched. 

It is worth to note that a different kind of process would 
lead to the same description. Consider for example a case 
where crystallite oxidation and corrosion are two possible 
processes in competition when a site is uncovered. If the 
probability p associated to the corrosion is proportional 
to the global etchant concentration, then a probability 
1 — p should be associated to oxidation (and then pas- 
sivation). In order to decide the etching or passivation 
of a given site i, a random number n is thrown. If 
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is smaller then p the site is corroded, otherwise is defin- 
tively passivated. The r<j numbers define a stochastic 
process which would give the same dynamical behavior. 
Of course the hardening properties would be different in 
that case. From a statistical point of view, it means sim- 
ply that we can equivalently formulate our model as a 
deterministic dynamics with quenched disorder or as a 
stochastic dynamics without quenched randomness. 



VII. CONCLUSIONS 

In this paper we have discussed several aspects of a 
simple model for the etching of a two-dimensional dis- 
ordered system by a finite volume of corroding solution. 
This has been done both theoretically and verified numer- 
ically . The dynamics correspond to the disappearance of 
weak surface sites which at the same time uncovers new 
sites. As the etching process consumes the etchant, the 
etching power of the solution decreases and the surface 
resistance increases to the point where the process stops 
spontaneously. One obtains a kind of "equilibrium" or 
static situation in which the dynamics is stopped. This 
static state is characterized by the fact that the surviv- 
ing interface sites have a resistance to the etching which 
is larger than the final value pf of the solution etching 
power which is on the order of the percolation threshold 

Pc- 

An analytical description of the time behavior of the 
solution etching power p(t) and of the distribution of re- 
sistances on the total interface has been introduced. This 
analytical approach indicates why and how the dynamics 
can be divided into two regimes. The first initial period 
corresponds to a classical, or superuniversal, regime. It 
can be described with precision by a mean field approx- 
imation. The second regime is a critical regime related 
to percolation criticality. The final connected interface 
is constituted by a collection of fractal interfaces up to 
a certain characteristic depth or scale a. The fractal di- 
mension is found to be very close to Df = 7/4. The 
difference p c — Pf between p c and the final etching power 
Pf, and the width a are both linked to the geometrical 
and external parameters characterizing the system via 
simple scaling relations. These properties can be simply 
explained by relating the model to the Gradient Percola- 
tion model: identifying the ratio L/V between the size of 
the solid and the volume of the solution with the gradient 
Vp which characterizes the scaling properties of GP. Af- 
ter this identification, it has been shown that our etching 
model belongs to the GP universality class, and that the 
exponents can be explained through percolation theory. 

An important result of this approach is the identifica- 
tion of the meaning of p(t) as the inverse of the mean 
number of new interface sites uncovered by each etched 
site. This identification in particular is very important 
in relation to the static final situation in which p(t) ~ p c , 
as it provides the physical meaning of the percolation 



critical threshold. 

Several further developments of these studies can be 
suggested. The statistics of other observable quantities, 
as the arrest time of the process tf or the maximal depth 
attained by the solution, can be studied and related to 
known results of asymptotic extreme theory Jl9| . In 
particular, such statistics determine the probability of 
"chemical" fracture of a finite solid submitted to an etch- 
ing process. Furthermore, the distribution of the debris 
produced by the etching process, can be regarded as a 
"chemical" fragmentation process. 

Also, the stability of the final (harder) interface with 
respect to external perturbations (as for example, fluctu- 
ation of the etching power pf) should be of interests. 
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APPENDIX A: THERMODYNAMIC LIMIT 

In light of what was written about the "phase" dia- 
gram, one can discuss the thermodynamic limit. Let us 
start with a couple of external parameters (Lq, Vq) in the 
strong gradient "phase" for which a(Lo, Vq) <C Lq, and 
negligible boundary effects. If one changes L (horizon- 
tal size of the solid) and V (volume of the solution) by 
satisfying the relation: 



L _ ( V \ 

the new system remains in the strong gradient "phase" 
for any value of V. The final corrosion front obeys the 
relation: 

o{L,V) = a{L ,V Q ) 
L L Q ' 

i.e. the new system is geometrically "similar" to the old 
one. 

Let us now study what happens changing (L, V) fol- 
lowing the relation 



L 

with a a a I (1 + a a ). 



(Al) 



• If a > 1 the system stays in the strong gradient 
phase for any value of V, but a(L, V) is a decreas- 
ing function of V and in the thermodynamic limit 
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(V — > oo), the final corrosion front is "microscopi- 
cally" flat; 

• If a a /(l + a a ) < a < 1 the system stays in the 
strong gradient phase. Moreover a(L, V) increases 
with V (it is infinite in the thermodynamic limit), 
but <j(L, V)/L < <t(L , Vo)/Lq. Hence the new cor- 
rosion front is not "similar" to the old one, and in 
the limit V — * oo it becomes "macroscopically" flat; 

• If a < a a /(l + a a ), a(L,V)/ L increases with 
V . However there will be a marginal (or "criti- 
cal") value V c of the volume for which the geo- 
metric correlations reach their maximum possible 
value: a(L, V) ~ L. Increasing V further, follow- 
ing Eq. [Al[ the system enters the weak gradient 
phase dominated by boundary effects. 



APPENDIX B: THE "GLOBAL" HISTOGRAM 

In the smooth regime, at each time step t, the number 
of surface sites which are created is L, and the number of 
dissolved sites is p(t)L. Consequently, at each time-step, 
the total number of surface sites increases by (1 —p(t))L. 
More precisely, we can see that these last sites affect only 
the high r part of the histogram. In fact, using Eq. ^5] in 
Eq. |2^, one can write the explicit form of h(r, t) for any 
time-step in the smooth regime: 



h(r, t)-- 
' L 



(Bl) 



for r < p(t — 1) 



= < (t - t')L for p(t' + 1) < r < p(t') ; t' < t - 2 



t [t + l)L for r > p(0) 



Note that h(r, t) is, at any time, a non-decreasing multi- 
step function of r. Each step differs from the previous 
one by height L. 

It can be written also: 



t-2 



h(r, t) = L6 (p(t - 1) - r) + ^ [(* - t')L> 



(B2) 



t'=0 



x9 (r - p(t' + 1)) 9 (p(t') - r)} + (t + l)L 9(r- p(0)) . 

In a similar way, the explicit function G(t) is obtained 
by inserting Eq. ^5| and Eq. ^ in Eq. ^2|: 



G(t) = L 



t + l-p(0) 



1 - cxp(-|) 
l-exp(-i) 



(B3) 



In order to obtai n th e nor maliz ed histogram c/>(r, t) one 
has to divide Eq 



shape of Eq. B2 



B2 by Eq. B3. Because of the step-like 
<fi(r, t) will also be step-like. 
A more physical derivation of a smooth function inter- 
polating the step-like </>(r, t) can be obtained under the 
same approximate phenomenological approach leading to 



Eq. 0. Under these assumptions, the corrosion front is lo- 
cated at depth y at the time t = y when the solution has 
the corrosive power p(t = y), where p{t) is given by Eq. || 
(or alternatively by Eq. ||(]) . From Eq. [| one can deduce 
that the solution attains the etching power p at the time 
t{p) given by 



*(?) = 



Hp/p(Q)) 

hx(l-L/V) 



(B4) 



when the front is at depth y = t{p). From this equation, 
it is possible to infer that a site with random resistance 
p(0) < r < p(t) is etched only if it is located at a depth 
y < t(r). In fact, if y > t(r), the site would be reached 
by the solution (i.e. checked by the dynamics) when the 
etching power p(t) is weaker than its resistance r. Hence, 
it is necessary to distinguish three cases, in order to write 
the number of solid sites, with resistance in [r, r + dr], 
belonging to the global surface at time t: 

1. All the sites with r > p(0) checked by the dynamics, 
resisted to the corrosion, and hence belong to the 
surface. Their number is given by dr multiplied 
by the area spanned by the front up to time t — 1 
(included), in addition to such sites on the corrosion 
front at time t: (Lt + L) dr; 

2. All the sites with p(t) < r < p(0) checked by the dy- 
namics before time t(r), have been etched, whereas 
such sites, checked between t(r) and t, resisted. 
The number is given by dr multiplied by the area 
spanned by the front between times t(r) and t — 1 
(included), in addition to the sites on the front at 
time t: [{L(t - t(r)) + L]dr; 

3. All the sites with r < p{t) checked by the dynamics, 
have been etched. Only sites with such resistence 
belonging to the corrosion front at time t contribute 
to the histogram. Their number is L dr. 

One can then write: 



h(r,t) 



L for r < p(t) 

L (t - t{r) + 1) for p(t) < r < p(0) (B5) 

L(t + 1) for r > p{0) 



Using the explicit formula for t(r) given by Eq. B4, with 
r replacing p, one has: 



G(t) = / h(r,t)dr = L 



t + 1- rp(0) 1 - e 



(0)(l 



(B6) 



where r = (^(l/L/U))" 1 ~ V/L. Note that Eq. [|| 
differs from the rigorous Eq. [B^ only from the approxi- 
mation exp(— 1/t) ~ 1 — 1/t which is valid in the present 
study where V/L 3> 1. The normalized distribution at 
all time is (f>(r,t) = h(r,t)/G(t), and for t ^> 1, it can be 
written as: 
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0M) 



01 (t) x 

for r < p(t) 
f (lnp(O) - lnr) forp(t) < r < p(0) 
for r > p(0) 




(B7) 



where 



[^(t)]- 1 - 1 + (p(0) -p(t))/ln(p(t)/p(0)) - 

= l-jp(0)(l-exp(-t/r)). (B8) 
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Eq. B7 is the same as Eq. ||l]. 
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